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ry, • Abstract. We study the Hubbard model at half band-filling on a Bethe lattice 

with infinite coordination number in the paramagnetic insulating phase at zero 
^5 ' temperature. We use the dynamical mean-field theory (DMFT) mapping to 

a single-impurity Anderson model with a bath whose properties have to be 
determined self-consistently. For a controlled and systematic implementation 
of the self-consistency scheme we use the fixed-energy (FE) approach to the 
DMFT. In FE-DMFT the onset and the width of the Hubbard bands are 
adjusted self-consistently but the energies of the bath levels are kept fixed 
w ^ relatively to both band edges during the calculation of self-consistent hybridization 

vJ , strengths between impurity and bath sites. Using the dynamical density-matrix 

renormalization group method (DDMRG) we calculate the density of states with 
a resolution ranging from 3% of the bare bandwidth W = 4t at high energies to 
0.5% in the vicinity of the gap. The DDMRG resolution and accuracy for the 
density of states and the gap is superior to those obtained with other numerical 
methods in previous DMFT investigations. We find that the Mott gap closes 
at a critical coupling Uc/t = A. 45 ± 0.05. At U = 4.5t, we observe prominent 

^ti* ' shoulders near the onset of the Hubbard bands. They are the remainders of the 

>0 ' quasi-particle resonance in the metallic phase which appears to split when the 

^_^ ' gap opens at Uc- 
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j^ ■ PACS numbers: 71.10.Fd, 71.27.+a, 71.30-fh 
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In the limit of higli dimensions Pj, models for correlated lattice electrons can be 
. , mapped onto effective single- impurity problems [3 E]. In some cases, the exact 

j^ ■ solution for a many-particle Hamiltonian has been found, e.g., for the Falicov-Kimball 

model 13 E], and for other problems a few exact results have been obtained; for 
reviews, see |S1 E]- Despite its increasing popularity [S] [7| |H1 it must be kept in 
mind that the dynamical mean-field theory (DMFT) still poses a difficult many-body 
problem: the effective single-impurity model must be solved self-consistently for the 
one-particle Green function at all frequencies. Consequently, reliable numerical or 
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analytical 'impurity solvers' must be developed and the self-consistency scheme must 
be implemented in a controlled way. 

The Hubbard model at half band-filling provides an ideal test case for the 
dynamical mean-field theory. It describes s-electrons with a purely local interaction of 
strength U and electron transfer matrix elements —t/vZ between Z —^ oo neighboring 
sites on a lattice. On the one hand, the model contains an interesting quantum phase 
transition between the paramagnetic metal and the paramagnetic (Mott-Hubbard) 
insulator [H] at a finite coupling. On the other hand, for a Bethe lattice with a 
semi-elliptic bare density of states of width W = At, perturbation theory to fourth 
order in U/W jJUj and to second order in W/U ^Jj have been carried out at zero 
temperature, against which approximate analytical and numerical techniques can 
be tested. In this way, merits and limitations of analytical methods (Hubbard-III 
approximation |12j . iterated perturbation theory |18| . local moment approach |14j ) 
and numerical techniques [exact diagonalization schemes (ED) J15l I16| . numerical 
renormalization group (NRG) [^] have been revealed, together with the difficulties 
in implementing the self-consistency scheme in numerical approaches. 

The latter problem results from the fact that numerical approaches work with a 
finite number of sites to represent the continuous bath coupled to the impurity site. 
Thus the energy resolution is necessarily limited by finite-size effects. Moreover, it is 
not clear a priori how one can define a self-consistency condition for the discretized 
impurity problem such that the self-consistent solution is approached in a smooth and 
controlled way in the limit of an infinite number of bath sites. [In other approaches 
to the impurity problem, such as quantum Monte Carlo (QMC) simulations, the bath 
is not discretized but the imaginary time has to be discretized, which leads to similar 
problems.] In a previous work ^1 ^^ this problem has been solved by the 'fixed- 
energy' approach to the dynamical mean-field theory (FE-DMFT): (i) a frequency 
interval / is split into subintervals le of equal length, whose mid-points ee give the 
energies of the bath sites, and the density of states is put to zero outside this interval /; 
(ii) the hybridization strengths between impurity and bath sites is determined self- 
consistently for these fixed energies eg. Within the fixed-energy approach to dynamical 
mean-field theory, the resolution of the frequency interval I improves systematically 
with system size rig, and an extrapolation ris — > cxd becomes meaningful. As has 
been shown in Refs. |l()l lll| . the FE-DMFT combined with exact diagonalization 
[FE-DMFT(ED)] with Ug < 15 provides a reasonable description of the metallic phase 
for U < OAW and of the Mott-Hubbard insulator for U/W > 1.2. 

With exact diagonalization finite-size effects are prominent in the interesting 
region of the metal-insulator transition, i.e., for U ~ W. Consequently, 
numerical approaches are required which overcome the limitation of the exact 
diagonalization technique. The dynamical density-matrix renormalization group 
method (DDMRG) ^HI treats large systems (here with up to Ug — 161 sites) very 
accurately. It is an extension of the standard density-matrix renormalization group 
(DMRG) |19[ 1201 to the calculation of dynamical correlation functions. For the 
computation of a continuous spectrum, DDMRG is more accurate than previous 
generalizations of the DMRG to dynamical quantities such as the Lanczos-DMRG |23 1 
or the correction- vector DMRG |25]. The DDMRG has been applied successfully to 
the single-impurity Anderson model [53| and to DMFT calculations for the metallic 
phase of the Hubbard model ^Hl ■ 

In this work, we present results for the Mott-Hubbard insulator on a Bethe 
lattice with Z -> oo neighbors obtained with FE-DMFT combined with DDMRG 
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[FE-DMFT(DDMRG)]. In Section H we specify the Hubbard model, the effective 
single-impurity Hamiltonian, the corresponding one-particle Green functions, and the 
self-consistency condition. We also recall the results from perturbation theory in 
1/U . In Section |31 we summarize important aspects of the fixed-energy approach to 
the dynamical mean-field theory, and the DDMRG impurity solver. Details can be 
found in Refs. ^JI^DES]- In Section 01 we display the density of states, the gap for 
single-particle excitations, the ground-state energy and the average double occupancy 
as a function of the interaction strength U in the Mott insulating phase found for 
U > Uc = 4.45(±0.05)t. A short summary and conclusions close our presentation. 
A sum rule for the ground-state energy of the single-impurity Anderson model at 
self-consistency is derived in [Appendix A| 

2. Definitions 

2.1. Hamiltonian 

We investigate spin-1/2 electrons on a lattice whose motion is described by 

where cj , c^ are creation and annihilation operators for electrons with spin a =t, | 
on site i. The matrix elements i" are the electron transfer amplitudes between sites 

i and j, and tj j = 0. Since we are interested in the Mott insulating phase, we consider 
exclusively a half-filled band where the number of electrons N equals the number of 
lattice sites L. 

For lattices with translational symmetry, t:^^ — t(i — j), the operator for the 
kinetic energy is diagonal in momentum space, 



where 



eik) = j-J2ti^-j)^''^'~'^'- (3) 

The density of states for non-interacting electrons is then given by 

k 

In the limit of infinite lattice dimensions and for translationally invariant systems 
without nesting, the Hubbard model is characterized by p{e) alone, i.e., higher-order 
correlation functions in momentum space factorize |24) . For our explicit calculations 
we shall later use the semi-circular density of states 



2 / f^LuV ,. . W^ 



W/2 

dujpo{uj) , (6) 

-w/2 
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where W = At is the bandwidth. In the foUowing, we take i = 1 as our unit of energy. 
This density of states is reahzed for non-interacting tight-binding electrons on a Bethe 
lattice of connectivity Z ^ oo [53] . Specifically, each site is connected to Z neighbors 
without generating closed loops, and the electron transfer is restricted to nearest- 
neighbors, %:? = —t/\/Z when i and j are nearest neighbors and zero otherwise. The 
limit .Z ^ oo is implicitly understood henceforth. 

The electrons are taken to interact only locally, and the Hubbard interaction reads 

^ - E (-U - \) {hi - 1) ' (7) 

i 

where n^ ^ — ct cj^ is the local density operator at site i for spin a. This leads us to 

the Hubbard model |26|. 

H = f + UD . (8) 

The Hamiltonian explicitly exhibits particle- hole symmetry, i.e., H is invariant under 
the particle-hole transformation 



C -(-l)l^lc. ; a. - (-l)l^lct , (9) 






where \i\ counts the number of nearest-neighbor steps from the origin of the Bcthe 
lattice to site i. The chemical potential /i = then guarantees a half-filled band for 
all temperatures |S]. 

2.2. Green Functions 

The time-dependent local single-particle Green function at zero temperature is given 
by IZa 

Git) = -^TJ2(nciAtHJ) ■ (10) 



■I, a 



Here T is the time-ordering operator and (. . .) implies the average over the degenerate 
ground states with energy Eq, and (taking h = 1 henceforth) 

c.Jt) = exp{iHt)c^^ eM-'^Ht) (11) 

is the annihilation operator in the Heisenberg picture. 

In the insulating phase we can readily identify the contributions from the lower 
(LHB) and upper (UHB) Hubbard bands to the Fourier transform of the local Green 
function {rj = 0+), 

/CXD 
dte'"*G(i) = Glhb(w) + Guhb(c^) , 
-CXJ 

Guhb(c^)= -Glhb(-^). (12) 

The last equality follows from the particle- hole symmetry Q . Therefore, it is sufficient 
to evaluate the local Green function for the lower Hubbard band which describes the 
dynamics of a hole inserted into the system. 
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The density of states for the lower Hubbard band can be obtained from the 
imaginary part of the Green function l|12|l for real arguments via j27j 

-Dlhb(^) = -QGlhbIw) , 



FEfe'5(- + ^-^")0' (13) 



L 

with uj < — A([/)/2 < 0, where A(C/) is the single-particle gap. Particle-hole symmetry 
results in a symmetric density of states around cu = at half band-filling 

D{u;) = i^LHB(^) + Dvhb{cj) (14) 

with DuhbI^) = -Dlhb(-w)- 

We define the (shifted) moments Af„([/) of the density of states in the lower 
Hubbard band via 

r-MU)/2 / rJ^n 

M„{U)=I du;iLJ + -j Di^HBiio). (15) 

In particular, from (|13|l we find that |2Z] 

Mo{U) = 1 , (16) 

1 / . X dEo(U)\ U , , 

Afi(C/) = -(^£;o([/) + C/^^j+- (17) 

are two useful sum-rules which we shall employ later. We also note that the average 
double occupancy is related to a derivative of the ground-state energy by 

2.3. Results from strong- coupling perturbation theory 

We shall test our numerical results against those from strong-coupling perturbation 
theory ^J. To second order in l/U, the density of states of the lower Hubbard band 
reads 

DLRBii^) = J ^ depo(e)s(e, U)S U + ^ + <?(e, U)\ + 0{U-') , 

e 9(e^-l) 
s{e,U) =1-^ + ^^, (19) 

. rn ^^'-3 3e(2e2-7) 

^^ ' ' 2U 8C/2 

The zeros of -Dlhb(w) provide the single-particle gap and the width of the Hubbard 
bands W*{U), 

A{U) =L/-4-l- J^+0(C/-3), (20) 

W*iU)=A+^+OiU-^), (21) 

up to second order in l/U. The Hubbard bands display a square-root onset, 

„„,,M.L_^)-" , „.^, ,22) 
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Note that there is no weight outside the Hubbard bands up to and including order 
1/U^ but there are contributions to order 1/U^ and higher. Our numerical results 
show that the weight outside the (primary) Hubbard bands at \llj\ > A{U)/2 + W*{U) 
is at most one percent of the total density of state for all interaction strengths in the 
insulating phase. 

Recently, E. Kalinowski [2111221 has calculated the ground-state energy to 11th 
order in the inverse coupling strength. Here we restate her results, 

EojU) ^_U_J- 1 ]^_^^_ 23877 _ 4496245 _ ^^(rj-is.foo. 

L 4 2U 2C/3 8t/5 32C/7 usu^ 2048C/ii ^ '^ ' 

1 3 95 4151 214893 49458695 ,^,^^ ,,, 

'^''^ = 2U^+2U-^ + m^+mi^ + T28mo+m8U^+''^''^^- ^''^ 
Unfortunately, the computational effort increases exponentially as a function of the 
order, and it will be difficult to obtain much higher orders of the expansion. 

3. Fixed-energy dynamical mean-field tiieory witii dynamical 
density-matrix renormalization group [FE-DMFT(DDMRG)] 

In this section, we first discuss the single-impurity model onto which the Hubbard 
model can be mapped in the limit of infinite dimensions. Next, we recall the 
fixed-energy algorithm for the dynamical mean-field theory. Lastly, we discuss the 
density-matrix renormalization group for the numerical solution of the single- impurity 
Anderson model. 

3.1. Dynamical Mean-Field Theory (DMFT) 

In the limit of infinite dimensions [J, and under the conditions of translational 
invariance and convergence of perturbation theory in strong and weak coupling, the 
Hubbard model can be mapped onto single-impurity models [H El Eli which need 
to be solved self-consistently. In general, these impurity models cannot be solved 
analytically. 

For an approximate numerical treatment various different implementations are 
conceivable. One realization is the single-impurity Anderson model in 'star geometry', 

-ffsiAM = X! ^i'^t-A'^;i + U ( d^d-i - 2 ) ( ^"t'^i ^ 2 

where Vi are real, positive hybridization matrix elements. The model describes the 
hybridization of an impurity site with Hubbard interaction to rig — 1 bath sites 
without interaction at energies eg. Here (i+, do-, ^^-^ , V'ct;£ are creation and annihilation 
operators for electrons with spin a =t,i on the impurity and the bath site t, 
respectively. In order to ensure particle- hole symmetry, we have to set e^ = — e„^_£ 
and Vi — Vn^-i for £ — {ug + l)/2, . . . , Tig — 1. Moreover, since we are interested in the 
Mott-Hubbard insulator, we only use odd Ug so that there is no bath state at e = 0. 
For a given set of parameters {eg, Vt] the model (|25|l defines a many-body problem 
for which the single-particle Green function 

GiSW--i(^[rf'.Wrf^])gj^^ (26) 
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must be calculated numerically. Here, (. . .)siam implies the ground-state expectation 
value within the single-impurity model. 

Ultimately, we are interested in the limit n^ -^ oo where 

Lo — ei + iO+ sgn(a;) 



^^"=H-) - E . .-,=^^.-....^ (27) 

i=i 

becomes the hybridization function of the continuous problem, 

H{lj) = lim H^"'\uj) (28) 

and the Green function is 

G,,iu;) = ^lim^[G(:r;(c.) + G^^^Jic)] . (29) 

[For finite Us the Green functions G^^£.^{u}) are different for a =],[ because the system 
contains an odd number n^ of electrons.] As shown in [2, the hybridization function 
and the Green function must obey a self-consistency relation, 

H{oj) = ^^^^ (30) 

to describe the Hubbard model on the Bethe lattice with connectivity Z ^ oo. At 
self-consistency, the Green function of the impurity problem gives the Green function 
of the Hubbard model, 

Gdd{i^) = G{lo) . (31) 

For a finite-size representation of the bath, n^ < oo, it is generally not possible to find 
a self-consistent solution to the finite-system version of H3l)|l . 

i?^""H'-) = ^[4:;:?H+G(':r](a;)]. (32) 

Instead, we have to choose bath energies e^ and hybridizations Vg for finite n^ in such a 
way that the single-particle Green function and the hybridization function fulfill (|30|l 
for ris — > 00. Therefore, numerical methods will differ in the way an approximate 
self-consistency condition is defined. This is a source of ambiguity because there can 
be more than one self-consistent set of parameters {ee, Ve} for fixed Us- Moreover, it 
cannot be guaranteed that different schemes will ultimately coincide for n^ — > cx). 

3.2. Fixed-energy dynamical mean-field theory (FE-DMFT) 

In Ref. m a new algorithm for solving the self-consistency problem has been 
introduced. The accuracy and stability of this 'fixed-energy DMFT' approach has 
been demonstrated using an exact diagonalization technique as 'impurity solver', i.e., 
to compute the single- impurity Green function GJ^''^{uj). In this work, we describe 
how to use the FE-DMFT together with the dynamical density-matrix renormalization 
group as impurity solver. 

For the Mott-Hubbard insulator, we make the assumption that all the spectral 
weight is concentrated in the upper and lower Hubbard bands, i.e., in the finite 
frequency interval 

I-^u. ^<H<^ + W*iU)\. (33) 
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The onset of the upper Hubbard band, A([/)/2, and the width of the Hubbard bands, 
W*{U), are determined self-consistently; see below. We start with some input guess 
A([/) and W*{U), which we may take from second-order perturbation theory H2Ufl . 
from the fixed-energy dynamical mean- field theory with exact diagonalization |11| , or 
from previous runs for slightly different values of U or n^. We discretize the Hubbard 
bands equidistantly, i.e., we fix the energies ee by 

ei^-en,-e = ^Y^+U-\]s{U) , 1 < £ < {n, - l)/2 , (34) 

where 

m) - =? (35) 

is the distance between two consecutive energies tt in the same Hubbard band. Then 
we divide the interval / into ris — 1 intervals Ii of width <5(C/) centered around each 
energy e^. By fixing the energies at the centers of equidistant intervals we can be sure 
that our resolution of the Hubbard bands becomes increasingly better as n^ increases. 
For a typical n^ = 65 and W*{\J) k. At we have d{U) « 0.125. 

When we integrate the imaginary part of the Green function over the interval Ig 
we obtain weights wi, 

Wi = / du . (36) 



At self-consistency ()30|l and for n^ ^ oo, these weights obey 

y/ = we . (37) 

We can use this relation to calculate new parameters Ve from a Green function Gddi^)- 
As initial values for Gdd{^) we may again use the result of second-order perturbation 
theory ^ in (EEJ, the results of the FE-DMFT(ED) [TTl, or we start from previous 
runs for slightly different values of U or n^. The latter approach is recommendable 
close to the transition. 

At every iteration the impurity Green function Gdd{^) nrust be calculated with the 
help of an 'impurity solver'. Here, we use the dynamical DMRG to calculate Gj^fllu;) 
for the Hamiltonian H25|l with finite Ug. Then, the deconvolution of the sum of these 
Green functions for a =1,1 gives an excellent approximation of the Green function 
Gdd{'^) in the limit n^ — > oo at all needed frequencies (see the next subsection). 

We now describe the iterative procedure used to determine the onset of the 
upper Hubbard band A(C/)/2, its width W*{U), and the Green functions Gdd{^) 
self-consistently. Starting from the initial A([/)/2, W*{U), and Gdd{^), we compute 
the energies ee and hybridization matrix elements Vi of the single-impurity Anderson 
model (|25ll using equations H34|) to (|37|l . In a first calculation we consider this model 
with Us sites and use the DDMRG method to compute the full Green functions 
G^ili^) with a resolution ry ~ <5(f/) ~ l/f^s- As explained above, after deconvolution 
of these Green functions we obtain a new Green function Gdd{^), which is used in the 
next iteration. Simultaneously, we use the DDMRG method with a broadening 77 <C 
6{U) to compute the energy A{U,n'^)/2 of the first pole in G^][l{uj) (i.e., the lowest 
state contributing to the density of states) for the single-impurity Anderson model 
with n'^ > Hs sites. Typically, we calculate A(C/, n'^) for n', = 81, 97, 113, 129, 145, 161. 
[These calculations can be carried out for larger system sizes than the calculation of 
the full Green functions because we only need to determine ground-state properties 
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and a small fraction of the Green function spectrum around lo ~ A(C/)/2.] After 
extrapolating to the thermodynamic limit, 

A(C/)= lim A(C/,0 , (38) 

we obtain a new estimate for the onset of the upper Hubbard band, A(t/)/2, which is 
used in the next iteration. At the same time we use the sum rule of [Appendix A| for 
the ground-state energy £^q^^^([/, n^) of the effective single- impurity Anderson model 
to calculate a new bandwidth, 

2 n^-^oo n'g 

After a new gap A([/), bandwidth W*{U), and Green function Gdd{'-^) have been 
obtained we can start the next iteration. We repeat this procedure until it converges to 
a fixed point. Typically, we need less than 10 iterations for the procedure to converge, 
depending on the choice of the starting parameters. We terminate the iterative 
procedure when the variation of the gap A(C/) and bandwidth W*{U) is smaller than 
lO^^i from one iteration to the next. At that point the variation of Gddi^) is found to 
be smaller than 10~'^ for all frequencies w. This iterative procedure is stable; for small 
deviations from the self-consistent values, the gap and the width of the Hubbard bands 
are driven back to the fixed point of the iteration. We have also checked that, for fixed 
Us, a unique solution for Gdd{^) is found for various starting choices. Obviously our 
FE-DMFT (DDMRG) approach yields self-consistent results for the gap, bandwidth, 
and Green function of the Hubbard model. Moreover, it is possible to calculate ground- 
state properties of the Hubbard model (energy, double occupancy) from ground-state 
properties of the self-consistent single-impurity Anderson model, as shown in the next 
section. 

In Figure n we give an example of the extrapolation scheme for A{U,n'^) at 
the fixed-point of our iterative procedure for U = 4.6. As expected, the results for 
81 < n'^ < 161 extrapolate linearly in l/n(.. Note that the FE-DMFT with exact 
diagonalization "11" was limited to n^ = n^ = 15, and finite-size effects had to be 
controlled by a combination with the criterion of a square-root onset of the Hubbard 
bands, which is suggested by perturbation theory 122|) . The DDMRG treats system 
sizes up to n'j, ~ 0(200) which makes this 'weight criterion' obsolete. Nevertheless, 
we can use the 'weight criterion' as a consistency check. As argued in Ref. |11) . we 
should find 

MU,0 f ^2/3 1,^. 

ei ^ ex (wi) , (40) 

for a square-root increase of the density of states near the band edges. In Fig. Q] we 
show ei as a function of Wi for system sizes 81 < n(, < 161 as open circles. Both 
extrapolated values for the gap from H38(l and (|40|l agree. The linear behavior of ei 
as a function of Wi confirms the square-root increase of the density of states near 
the gap. Note, however, that the region in which the square-root onset is discernible 
becomes very small close to the transition and thus large system sizes are required. 

For U < 4.6, a constant discretization of the Hubbard band with S{U) « 0.125 is 
not sufficient to resolve fine structures of the density of states near the single-particle 
gap. In order to obtain a better resolution for |w| « A([/)/2 without excessive increase 
of the computational effort we use a variable discretization scheme as described in 
Ref. |23|. The resolution around the gap is improved by using a finer discretization 
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Figure 1. Extrapolation of the lowest lying single-particle excitation energy 
A{U,n'g)/2 as a function of inverse system size (lower axis) for U = 4.6 at self- 
consistency (solid circles). The open circles give the energy ei as a function of the 
weight wi (upper axis) for the same system sizes n'^ = 81, 97, 113, 129, 145, 161. 



6{U) of the intervals A{U)/2 < \u!\ < t (i.e., more bath states are used in those 
intervals). The smaller 6{U) allows us to use a smaller broadening r/ in DDMRG 
calculations for those frequencies. We combine the high-energy spectrum obtained 
with the usual resolution and the low-energy spectrum obtained with the improved 
resolution and then deconvolve the result to obtain a new Green function Gdd{^)- This 
yields Gddi'^) for |cl;| < 0.6t with a resolution, which is up to an order of magnitude 
better than for a constant discretization with ris = 65. For the results presented here 
we have used 6{U) = 0.02 (corresponding to Ug = 113 and 77 = 0.03) for U = 4.5 and 
d{U) = 0.031 (n^ = 97 and rj = 0.05) for U = 4.6 in the intervals A{U)/2 < \uj\ < t. 



3.3. DDMRG for the single-impurity Anderson model 

The DDMRG for the single-impurity Anderson model is described in detail in Ref. [23 ■ 
Here, we summarize the essential ingredients. As the DMRG method is most accurate 
for systems with a quasi one-dimensional structure, we perform calculations of the 
single-impurity Anderson model 125() in its equivalent linear-chain form j3()| 

i^siAM - U (d+d^ - I) (d+d^ -l)+vy^ (f+J, + d+f,,o) (41) 
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a £=0 

The DDMRG provides the local density of states 



o'LM) 



-sgn(wi)- 



(42) 



at selected frequencies uji very accurately. The real part of the Green function can 
also be calculated with DDMRG but to carry out the FE-DMFT calculation we need 
only the imaginary part. To simulate the continuous spectrum of an infinite chain 
in a calculation with a finite n^, a broadening r\ is introduced which must be scaled 
as a function of the system size JHl- If *? is chosen too small, the density of states 
displays finite-size peaks as those seen in Ref. |3J|- If ^7 is chosen too large, relevant 
information is smeared out. As an empirical fact, r\ ~ <5(t/) = 2W^*([/)/(ns — 1) 
should be chosen, i.e., the resolution scales as the inverse system size, as found for 
one-dimensional lattice models I18j . 



gO.5 




u 



4.6 4.7 4.8 4.9 



Figure 2. Main figure: density of states of the upper Hubbard band for C/ = 6. 
The dashed Une shows the result of the DDMRG with a broadening rj = 0.2. 
The circles give the DDMRG density of states after deconvolution (77 = 0). The 
full line is the result from second-order perturbation theory in 1/U I19i . The side 
figures show the linear behavior of the square of the density of states as a function 
of frequency near the band edges. The lines are linear fits. 



In order to carry out the iterative procedure described in the previous section, we 
determine the density of states at selected frequencies w^. Typically, we choose them 
to resolve the effective bandwidth W*{U) equidistantly, uji^i — uji = Suj ~ rj ^ 5{U). 
We then 'deconvolve' the DDMRG data by inverting the Lorentz transformation [TUI 

5lo 



DU^i) - 



V 



j 



TT if + [uJi - LUj)' 



;Ddd{uJj) 



(43) 



where D^^{u!) — D^^ Ju!) + D^^ li^)- Through equation (|42|l this deconvolved density 
of states Dddifjj) determines the imaginary part of the Green function Gdd{^) which 
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is used in the FE-DMFT (DDMRG) scheme. The procedure can be repeated for 
different choices of the equidistant frequencies loj to get more values of D4^,j{u)j). 
In practice, we use two different sets of frequencies, corresponding to a frequency 
resolution comparable to the bath energy resolution 5{U). In this way, DDMRG 
provides a set of values Ddd,a{^j) for the density of states. The main advantage 
of this deconvolution is that no extrapolation or finite-size scaling analysis of these 
values Ddd,a{uJj) is necessary because they converge very quickly to the ns ^ oo limit. 
Naturally, structures with an intrinsic width of less than 77 cannot be resolved with this 
procedure even if we use many different sets of frequencies. Therefore, with DDMRG 
we obtain an accurate discrete representation of the density of states for the continuum 
model [and thus of the imaginary part of G'dd(w)], except for small regions around its 
onset and closing points where the derivative of the density of states changes singularly. 
With the DDMRG method ^HIj we calculate the one-particle Green function (|26|l 
for system sizes up to n^ ~ 0(200). Therefore the FE-DMFT (DDMRG) leads to a 
much better resolution of the Hubbard bands than our previous FE-DMFT with exact 
diagonalization which was limited to n^ — 15. 

An example of the density of states obtained with the FE-DMFT(DDMRG) 
approach is shown in FigureElfor C/ = 6. For this interaction strength, the agreement 
of the deconvolved DDMRG data with the second-order strong-coupling perturbation 
theory H19|l for the Hubbard model is almost perfect. Our deconvolution scheme 
gives slightly negative values in the vicinity of the band edges. These effects are 
small and are to be expected for sharp band edges in the density of states at 
w = A(C/)/2 and uj = A(C/)/2 + W*{u). We note that our numerical results 
are in much better agreement with perturbation theory than the results obtained 
in a recent DMFT(DMRG) study [HI where Lanczos-DMRG and a different self- 
consistency scheme has been used. Therefore, we think that our results for the gap 
and the critical interaction strength are also more accurate than those presented in 
that work pTTj . 

4. Results 

In this section we present the results for the Mott insulating phase of the Hubbard 
model which we have obtained with our FE-DMFT (DDMRG) approach. For ground- 
state properties comparisons with strong-coupling perturbation theory ^| 1231221 and 
DMFT(QMC) results (extrapolated to zero temperature) |28IIH2| confirm the accuracy 
of our method. Moreover, we will present results for the (zero-temperature) single- 
particle excitations which are much more accurate than those obtained with other 
DMFT approaches. 

^.1. Ground- state properties 

The ground-state energy per site of the Hubbard model can be calculated from 
ground-state expectation values of the self-consistent single-impurity Anderson model 
(see [Appendix A) ) 

^^ + J = U{d+d^d+d{)siAM + V{dtka)siAM , (44) 

where the two terms on the right-hand-side are the interaction and kinetic energy 
per site, respectively, and V = t = 1. In Figure 13 we show the ground-state energy 
Eo(U)/L + U /A in the Mott-Hubbard insulator phase for 4.5 < [/ < 6 in comparison 
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with strong-coupling perturbation theory H23|l . We see that there is a very good 
agreement between our numerical data and the analytical results. Our data points lie 
below the best perturbative energy (llth-order in l/U). As expected, deviations from 
the perturbative results become larger when U becomes smaller, from about 0.8 x 10"^ 
at t/ = 6 to 8.8 X 10-"* at f7 = 4.5. Our FE-DMFT(DDMRG) energies are also lower 
than DMFT(QMC) energies |28l IH2j . However, the differences are small, of the order 
of 2 X 10"'* or less, for U > 4.8. As the Mott insulator solution disappears for U < 4.8 
in the DMFT(QMC) approach, no comparison with our data is possible below that 
coupling strength. 



-0.08 




1^ order 
3^^ order 
7 order 
1 1 order 
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Figure 3. Ground-state energy Eq/L + U/4 of the Mott-Hubbard insulator 
as a function of the interaction strength. FE-DMFT(DDMRG) results for 
U = 4.5, 4.6, 4.7, 4.8, 5, 6 (circle) and perturbation theory (lines) for various orders 
in l/U. 



It is difficult to evaluate the relative accuracy of our method from the ground- 
state energy alone because that quantity is only defined up to a constant. The average 
double occupancy of the Hubbard model is given by the average double occupancy of 
the impurity site in the single impurity Anderson model at self-consistency 

d{U) = {d+d^d+ di)siAM . (45) 

At half filling this quantity takes only values between zero and 1/4 and thus provides 
a good benchmark. In Figure 01 we compare our FE-DMFT (DDMRG) results for 
the average double occupancy with perturbation theory (|24|l up to 12th order in l/U. 
Again, the agreement is very good but deviations become clearly noticeable for U < 5. 
Quantitatively, the differences between our values for d{U) and the results of the 12th- 
order perturbation expansion increase significantly from 2 x 10^^ (about 0.01%) at 
C/ = 6 to 2 X 10"'^ (about 7%) at U = 4.5. This is not surprising because we locate the 
critical value below which the Mott insulator no longer exists at Uc « 4.45 ± 0.05 (see 
below). The series expansion for the ground-state energy (^5)) and the average double 
occupancy (|24|l converges only for U > Uc. Therefore, the results of finite-order 
perturbation theory rapidly become inaccurate as U approaches Uc. A comparison 
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Figure 4. Average double occupancy in the Mott— Hubbard insulator as a 
function of the interaction strength U. FE-DMFT(DDMRG) results for U = 
4.5, 4.6, 4.7, 4.8, 5, 6 (circle) and perturbation theory for various orders in 1/f/ 
(lines). Inset: same results as a function of {U — U{)^''^ with C/f = 4.419. The 
dotted line is a fit to 1461 . 



between FE-DMFT(DDMRG) and DMFT(QMC) |2HI OH data is more conclusive. 
Both approaches provide results for the average double occupancy which deviate from 
each other by less than 3 x 10~^, corresponding to relative errors of the order of 0.1%, 
down to [/ = 4.8. 

With our FE-DMFT (DDMRG) approach the Mott insulator is stable for 
significantly weaker couplings than predicted by other works ^| EHl- A closer 
inspection of our data for small f7 < 5 shows that the double occupancy scales as 

d{U) =df- C^JU - Uf . (46) 

This behavior is clearly seen in the inset of Fig. ^ A fit of our data for [/ < 5 
gives Ui = 4.419, df = 0.03931 and C = 0.0202. Equation ^ suggests that the 
double occupancy is a singular function of the coupling [/ at C/f . It is thus reasonable 
to identify U[ with the critical coupling below which the Mott insulator no longer 
exists. The value U{ — 4.419 is indeed in very good agreement with the coupling 
Uc = 4.45 ± 0.05 where the Mott gap A(C/) closes (see below). As the average double 
occupancy is related to the ground-state energy by (|18|l . one expects that 

-^ + - = eo + d,{U -Ui)- — 

for U ^ Ui. Our data for the ground-state energy for [/ < 5 are reproduced by 
this formula within 5 x 10^^ if we use eo = 0.12235 and the parameters di^Ui 
and C determined from the fit of the double occupancy data. Therefore, our FE- 
DMFT(DDMRG) data for the ground-state energy and the average double occupancy 



-iU-Ui)^ 



(47) 
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of the Hubbard model with 4.5 < f/ < 5 fulfiU the relation (|18|l very precisely. For 
an arbitrary single-impurity Anderson model the derivative of the expectation value 
in the right-hand-side of H44|l is not given by the average double occupancy l|45|l . 
The relation H18|l between ground-state energy and double occupancy is valid for the 
Hubbard model and thus only for the expectation values H44|l and H45|l of single- 
impurity Anderson model at self- consistency. Therefore, the scalings H46(l and (|47|l 
of our data confirms that we have found self-consistent DMFT solutions 1)30(1 for the 
Hubbard model with couplings 4.5 < [/ < 5. 

Moreover, these results show that our FE-DMFT(DDMRG) approach is accurate 
enough to allow for an analysis of the critical behavior and to determine critical 
exponents for the ground-state energy and the average double occupancy. Note that 
the behavior H46(l of the average double occupancy implies that the interaction energy 
Ud{U) scales as C/fdf — CUi^/U — U[ close to U[ and, consequently, the kinetic energy 
scales as K+CUi\/U — U{ for small U—U{, where eo = K+Ufdf. Recently, evidence for 
half-integer critical exponents have also been found using an analysis of the strong- 
coupling perturbation theory extrapolated to infinite order |28| . However, the first 
singular terms in the expansions of Eo{U) and d{U) were found to have exponents 5/2 
and 3/2, respectively, compared to 3/2 and 1/2 in H47|l and (|46|l . 







4.4 4.6 4.8 
U 
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Figure 5. Single-particle gap in the Mott-Hubbard insulator as a function 
of the interaction strength. FE-DMFT(DDMRG) results (circle), second-order 
perturbation theory (solid line), and the interpolated result from FE-DMFT(ED) 
(dashed line). 



4-2. Single-particle excitations 

The single-particle gap A(C/) found at the fixed-point of our iterative procedure is 
shown in Figure |S1 As expected, A(C/) first decreases monotonically with U then 
vanishes below a finite Uc > 0. For U = 4.5 the gap A{U) = 0.062 is still large enough 
to be detected with our method but for U = 4.4 we find A(C/) = 0. Thus we estimate 
that Uc ~ 4.45 with an error smaller than 0.05 in full agreement with the singular 
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behavior of the ground-state energy (|17|) and double occupancy (jl^ described above. 
In Fig. it is seen that second-order perturbation theory describes the gap behavior 
quahtatively but it predicts a vanishing of the gap at a shghtly too smaU [4 = 4.31. 
We also see that our FE-DMFT(DDMRG) results agree with the resuhs from the FE- 
DMFT(ED) investigation ^I]. The small deviations are within the error estimates for 
FE-DMFT(ED) calculations (see Ref. 52)- Concomitantly, the values for the closing 
of the gap are almost equal, Uc = 4.43 ± 0.05 with the FE-DMFT(ED) method. 

Our result for Uc is in conflict with the value Uc = 4.78 found using a 
DMFT(NRG) approach ^7] and an analysis of the strong-coupling expansion '^. In 
contrast, we find substantial gaps A(C/ = 4.8) = 0.356 and A(f7 = 4.7) = 0.260 just 
above and below that coupling. These gaps are clearly larger than the discretization of 
the bath 5{U) = 0.125 that we have used. Therefore, we are confident that Uc < 4.7, 
and that our result Uc « 4.45 is more reliable than the results of Refs. |17ll28j . 

For large interaction strengths, the derivative of the gap A([/) with respect to U 
is unity, A'(C/ > VF) = 1, see (EOJ. For finite U, our results show that A'(C/) > 1, in 
agreement with perturbation theory H20() . In the vicinity of the transition, U ~ Uc, 
A'(t/) again approaches unity and thus A(J7) = U —Uc- 

The width of the Hubbard bands W*{U) calculated at the fixed-point of the FE- 
DMFT(DDMRG) procedure is almost constant for all U > Uc- At finite coupling 
it is slightly larger than the value W*{U) — W — A predicted by strong-coupling 
perturbation theory ioi U -^ 00 and reaches a maximum W*{U) « 4.1 around U = 5. 




Figure 6. Density of states of the upper Hubbard band for 4.6 < U < 5. 



In order to explain this behavior, we display the density of states as a function 
of U in Fig. |H| For large interaction strengths, U > 6, second-order perturbation 
theory describes the density of states D{uj) accurately, as seen in Fig. [3 The 
spectrum consists of the two Hubbard bands around ±f7/2 with a square-root onset 
at w = ±A(C/)/2. For weaker couphngs our FE-DMFT(DDMRG) calculations show 
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clearly that a shoulder forms in the density of states near the transition to the metallic 
phase. In Fig. we can see that this feature becomes progressively stronger as 
U approaches Uc- Its appearance is connected with the non- monotonous behavior 
of W*{U) and of A'([/) as a function of U near U — 5. 

This feature is the remainder of the quasi-particle peak which is present a.t ui — 
in the metallic phase. Because the metal is a Fermi liquid, the quasi-particle peak 
has height D{cj = 0) = p(0) = I/tt [^SI- Figure suggests that the quasi-particle 
peak splits at the transition to the insulating state at Uc- (A splitting of the density 
of states at the transition also occurs in the one-dimensional Hubbard model where 
Uc = 0+, as shown in Ref. ^J within a field-theoretical approach.) As the gap opens, 
the two flanks of the peak quickly loose weight so that they are rather small already 
at [/ = 4.5. 




Figure 7. Density of states of the lower and upper Hubbard bands for U = 4.5. 



Clearly, our FE-DMFT(DDMRG) results for the gap and the bandwidth are 
more accurate and the density of states shown in Figs. El to have a much better 
resolution than those calculated with other methods such as the FE-DMFT(ED) |llj . 
DMFT(NRG) ^7\ or DMFT(DMRG) EI]; DMFT(QMC) calculations 28 have not 
provided estimates for these quantities. In particular, our investigation demonstrates 
the presence of a sharp feature just above the gap in the density of states of the 
insulator and, thus, provides the first clear evidence for a splitting of the quasi-particle 
peak at the metal-insulator transition. As an accurate description of the low-energy 
excitations is necessary close to the critical coupling Uc, it is not surprising that 
the insulating phase extends to slightly weaker couplings than predicted in previous 
works ^1 1281 m] . The impurity solver used in previous DMFT investigations do not 
have the accuracy of the DDMRG method combined with the 'fixed-energy' dynamical 
mean-field theory. Therefore, they could not resolve the small gap A([/) < 0.260 for 
U < 4.7 and did not find a stable insulating phase below that coupling. 



5. Summary and conclusions 

We have investigated the paramagnetic insulating ground state of the Hubbard model 
on a Bethe lattice in the limit of infinite coordination number. In this limit, the 
problem can be treated within the dynamical mean-field theory (DMFT), i.e., it 
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can be mapped onto a system made of a single impurity with Hubbard interaction 
and hybridizations to a bath. The system properties have to be determined self- 
consistently from the required equivalence between the single-particle Green function 
and the hybridization function. In this work, we have used the fixed-energy approach 
to the DMFT (FE-DMFT) [TTj. The FE-DMFT provides stable solutions of the 
DMFT self-consistency problem and a systematic convergence of the results with 
increasing system size. 

As 'impurity solver' for the single-particle density of states we have used the 
dynamical density- matrix renormalization group (DDMRG) ^3- Our results from 
FE-DMFT (DDMRG) for the ground state agree with perturbation theory in l/U, 
where the latter is applicable, and with quantum Monte Carlo (QMC) simulations, 
where the DMFT(QMC) approach finds a stable insulating phase. With DDMRG we 
have used up to 160 sites to represent the self-consistent DMFT bath as compared to 
Hs = 15 with exact diagonalization (ED). These larger systems provide an enhanced 
resolution which is necessary to reveal structures in the density of states near the gap. 
These structures emerge when the Mott gap becomes small, i.e., when the coupling 
U approaches a critical value Uc below which there is no insulating phase (= Uc^i in 
a scenario with coexisting metallic and insulating phases). Our FE-DMFT(DDMRG) 
study gives Uc = 4.45 ± 0.05 for the critical interaction strength where the gap closes, 
in very good agreement with our previous FE-DMFT(ED) study JT], Uc = 4.43±0.05. 

In contrast to the results of a recent DMFT(DMRG) work |3J|, our results are not 
dominated by finite-size effects. At [/ = 6, for example, the density of states in [HJ 
displays a series of individual peaks instead of the smooth Hubbard bands found in 
our approach and in perturbation theory. Preliminary results for the metallic Fermi- 
liquid phase just below Uc suggest that the narrow quasi-particle resonance simply 
splits at Uc- Narrow shoulders which can be seen in the insulator density of states for 
U = 4.5 in Fig. [T] are the remainders of the quasi-particle resonance. The shoulders 
quickly loose weight as the upper and lower Hubbard bands separate from each other 
with increasing interaction strength U. 

The method presented here can also be applied to the metallic phase, as done 
in Ref. ^Uj for the weak-coupling limit. It is more difficult to resolve sharp quasi- 
particle peaks with DDMRG [22| than with, e.g., the numerical renormalization 
group. However, as shown in this work, our method offers the unique advantage 
that we can resolve sharp structures in the vicinity of the Hubbard band onsets. 
This is very important to describe the Mott insulating phase accurately and to 
determine the parameter region where it exists. Therefore, we are confident that 
our FE-DMFT (DDMRG) will provide deeper insight into the Mott-Hubbard metal- 
to-insulator transition. 
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Appendix A. Sum rule 

At self-consistency we liave from equations H13|l . (I15II . (|16|l . and (|3()|l 

Mi([/) 

cr '^LHB 

HslAM-Ef'^^)d^ (Al) 

/ /SIAM 



where £'q is the ground-state energy of the single-impurity Anderson model. The 
fact that the first moments are identical at self-consistency also proves that the 
average double occupancy of the Hubbard model d{U) is identical to the average 
double occupancy of the interacting site in the impurity model. Therefore, we do not 
distinguish between d{U) and (isiAM(f^)- 

We carry out the commutators in (jA.ip using the Hamiltonian (|25f) and obtain 

Afi([/) = 2Ud{U) -f Y. yddti^t,a)siAM . (A.2) 

The ground-state expectation value on the right-hand-side of this equation is readily 
calculated in DMRG. Therefore, combining (|17|l . I|18|) and l|A.2|l provides the ground- 
state energy density of the Hubbard model in terms of the single-impurity results 
as 

M^ = Ud{U) - J + E yi{dti^i,'rhiAM , (A.3) 

which is equivalent to (|44|l . 

The static ground-state expectation value in IjA.ip can be obtained from the 
corresponding Green function. For completeness we define the time-ordered Green 
functions 



Gu-At)^ -i(r[V'£..(t)7A+j) 



SIAM 



Gu-At) = - KT'[i^tAt)dt])siAM , (A.4) 

Gdi;ait) = -i(T[(io-(i)V'£'^cr])siAM • 

With the help of the equation of motion it is not difficult to show that their Fourier 
transforms obey {uj = ui + iO+sgncj) 

LU - €£ (u) - eg)'^ 



Gu-a{A — ~ ~ + 7~ —y^Gdd-cr{^) , 



Vg 

Ged:a{A = Gdd-a{oj) = Gde-a{A ■ (A. 5) 

uj — eg 

Then, the first moment in ljA.21) becomes {v = 0+) 

/OO J Tr 

^e'^-^^GddA^) ■ (A.6) 

At self-consistency we have (|30ll which implies 

E^H = E-^ = GH- (A.7) 
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Therefore, we arrive at the important relation 

^e'"" [G{uj)f = 2 (Afi([/) - 2Ud{U)) . (A.8) 

It relates the ground-state energy of the single-impurity model back to the ground- 
state energy of the Hubbard model. To show this we start from 

+ y^ ee{'ipl^il^e,cr)siAM ■ (A. 9) 

With the help of ljA.5|) and (|A.7p we find for the second term in (|A.9|) 

E ^^("^tJ- + V'^,-^^a >siAM = 2 (Afi((7) - 2C/d(;7)) . (A.IO) 



For the third term we use 






to obtain 



Y. ''Sta^e,a)siAM = ^ '^('^(-^e) (A.12) 



t,a 



2.1 ^ 2m ^ ' 






= ^ Qe(-e,) - i (Mi(C/) - 2Ud{U)) . 

Summing the contributions from l|A.10(l and l)A.12|l in 1A.9I) gives the final result 

Ef^^{U) ^ ^Qe(-Q) - 2Ud{U) + \m,{U) - ^ , (A.13) 

where d{U) is the average double occupancy l|18|l . Mi{U) is the first moment (|17f) . 
and 0(x) is the step function. This equation expresses the fact that the impurity 
provides corrections of order unity to the extensive ground-state energy. Therefore, 
for our equidistant energy levels eg, we find equation H39|). 
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